c
c     Second Derivative of a function
c     ------------------------------
      subroutine deriv_2 (f, n, d, h, IY)
      real*8 f(n), d(n), h, h2
c     f(n) - function
c     d(n) - It's second derivative
c     h    - step by x
c     IY .eq. 1 then use  error=h**2 - formular
c     IY .ne. 1 then use  error=h**6 - formular
      integer i, n, IY
      if (IY .eq. 1) then
         h2 = 1d0 / (h * h)
         do i = 2, n - 1
            d(i) = (f(i-1) + f(i+1) - 2d0 * f(i)) * h2
         end do
      else
         h2 = 1d0 / (12d0 * h * h)
         i = 2
         d(i) = (11d0 * f(i-1) - 20d0 * f(i) + 6d0 * f(i+1) +
     >      4d0 * f(i+2) - f(i+3)) * h2
         i = n - 1
         d(i) = (11d0 * f(i+1) - 20d0 * f(i) + 6d0 * f(i-1) +
     >      4d0 * f(i-2) - f(i-3)) * h2
         do i = 3, n - 2
            d(i) = (-f(i-2) + 16d0 * f(i-1) + 16d0 * f(i+1)  - f(i+2)
     >         - 30d0 * f(i)) * h2
         end do
      end if 
c
      h2 = 1d0 / (12d0 * h * h)
      i = 1
      d(i) = (35d0 * f(i) - 104d0 * f(i+1) + 114d0 * f(i+2) -
     >   56d0 * f(i+3) + 11d0 * f(i+4)) * h2
      i = n
      d(i) = (35d0 * f(i) - 104d0 * f(i-1) + 114d0 * f(i-2) -
     >   56d0 * f(i-3) + 11d0 * f(i-4)) * h2
c
      return
      end
c
c
c     First Derivative of a function
c     ------------------------------
      subroutine deriv_1 (f, n, d, h, IY)
      real*8 f(n), d(n), h, h2
c     f(n) - function
c     d(n) - It's first derivative
c     h    - step by x
c     IY .eq. 1 then use  error=h**2 - formular
c     IY .ne. 1 then use  error=h**4 - formular
      integer i, n, IY
      if (IY .eq. 1) then
         h2 = 1d0 / (2d0 * h)
         do i = 2, n - 1
            d(i) = (-f(i-1) + f(i+1)) * h2
         end do
      else
         h2 = 1d0 / (12d0 * h)
         i = 2
         d(i) =  (-3d0 * f(i-1) - 10d0 * f(i) + 18d0 * f(i+1) -
     >      6d0 * f(i+2) + f(i+3)) * h2
         i = n - 1
         d(i) = -(-3d0 * f(i+1) - 10d0 * f(i) + 18d0 * f(i-1) -
     >      6d0 * f(i-2) + f(i-3)) * h2
         do i = 3, n - 2
            d(i) = (f(i-2) - f(i+2) - 8d0 * (f(i-1) - f(i+1)) ) * h2
         end do
      end if 
      h2 = 1d0 / (2d0 * h)
      i = 1
      d(i) =  (-3d0 * f(i) + 4d0 * f(i+1) - f(i+2)) * h2
      i = n
      d(i) = -(-3d0 * f(i) + 4d0 * f(i-1) - f(i-2)) * h2
      return
      end
